Phonon effects in molecular transistors: Quantal and classical treatment 
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We present a comprehensive theoretical treatment of the effect of electron-phonon interactions 
on molecular transistors, including both quantal and classical limits. We study both equilibrated 
and out of equilibrium phonons. We present detailed results for conductance, noise and phonon 
distribution in two regimes. One involves temperatures large as compared to the rate of electronic 
transitions on and off the dot; in this limit our approach yields classical rate equations, which are 
solved numerically for a wide range of parameters. The other regime is that of low temperatures 
and weak electron-phonon coupling where a perturbative approximation in the Keldysh formulation 
can be applied. The interplay between the phonon-induced renormalization of the density of states 
on the quantum dot and the phonon-induced renormalization of the dot-lead coupling is found to 
be important. Whether or not the phonons are able to equilibrate in a time rapid compared to the 
transit time of an electron through the dot is found to affect the conductance. Observable signatures 
of phonon equilibration are presented. We also discuss the nature of the low-T to high-T crossover. 



I. INTRODUCTION 

In recent years it has become possible to fabricate 
devices in which the active element is a very small or- 
ganic moleculei. Such a device may be thought of as a 
'quantum dot': a structure weakly coupled to the macro- 
scopic charge reservoirs ('leads') and small enough that 
the discrete nature of the energy levels on the dot is 
important. Quantum dots fabricated using conventional 
semiconductor technology have been extensively studied 
experimentally 2 and theoretically 3 *. However, the use of 
small molecules may lead to new physics. In particular, 
as electrons are added or removed from a small molecule, 
both the shape of the molecule and its position relative 
to the leads may be altered. The energies associated with 
these changes are not small, and the time scales may be 
comparable to those related to the flow of electrons into 
and out of the molecule. Interesting recent data indicate 
that these effects may lead to observable structures in 
the conductance spectra of the dol££&. 

The shape change may be thought of as a cou- 
pling of electrons on the molecule to phonon modes of 
the molecule, while the position change corresponds to 
phonon dependent tunneling matrix elements. The sub- 
ject of electron-phonon coupling in quantum dots has 
received much theoretical attention 7 ! 8 ! 9 ! 10 ! 11 ! 12 ! 13 ! 14 ! 15 ! 16 . 
In an early important work Glazman and Sh.ckh.ter ob- 
tained analytic expressions for the transmission probabil- 
ity through the dot under conditions that the phonons are 
always in equilibrium 7 . Their results for the transmis- 
sion go a long way towards describing the behavior of the 
phonon-coupled system when one is far from resonance. 
However their treatment neglects the phonon renormal- 
ization of the dot-lead coupling, and thus gives rise to a 
zero bias conductance at resonance that is smaller than 
the value predicted by the Breit-Wigner formula (this 
situation was not the main interest of 7 ). (These issues 
were also recently discussed by Flensbergi 5 .). The lack of 
renormalization of the dot-lead coupling appears in the 



treatment carried out by other authors as well 8 *2ii2ii4 ) 
and in addition some authors assert (incorrectly, we be- 
lieve) that phonon sidebands may be observed even in 
the linear-response conductanc e) 9 ' 10 by tuning the gate 
voltage. 

In this paper we revisit the problem of the phonon cou- 
pled dot. We present a comprehensive formalism valid 
both in the classical and quantal limits which resolves 
the ambiguities in the present literature. We also use 
this formalism to address new issues related to the be- 
havior of this system under strongly non-equilibrium con- 
ditions in the quantal and classical regime. We further 
study the nature of the quantal-classical crossover and 
the noise spectra. An important feature of our results in 
the quantum regime is that the conductance peak height 
at resonance is unchanged by electron-phonon interac- 
tions. 

A recent paper by McCarthy et aliii treats the problem 
of the phonon coupled dot in the high-T classical regime 
and for phonons that couple to the leads. They present 
results for the regime of "equilibrated" phonons strongly 
coupled to a heat bath. The high-T regime of our work 
is similar to that of reference^, but we also study the 
physics of out of equilibrium phonons. 

Aji et al>iS have studied the model under off-resonant 
conditions when the conductance is very low and the cur- 
rent is due to elastic/inelastic cotunncling. They study 
the phonon sidebands in the case of equilibrated and un- 
cquilibrated phonons, for the case of an electron coupled 
to a molecular vibrational mode (our model) and also 
the case of a phonon dependent tunneling amplitude. In 
our analysis, the exact treatment of the leads automati- 
cally takes into account cotunneling processes, while the 
perturbative approximation restricts our analysis to de- 
tails of the first phonon sideband. We show how thermal 
effects wash out cotunneling. 

The paper is organized as follows. In Section II we 
describe the model (sub-section A), the important model 
parameters (sub-section B) and develop a density matrix 
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formalism (sub-section C) that allows one to obtain the 
probability distribution for various states on the dot un- 
der out of equilibrium conditions. In Section-Ill we ap- 
ply the high temperature approximation in the density 
matrix formalism and derive the rate equations for the 
dot occupation probabilities in the sequential tunneling 
regime (sub-section A). We use these rate equations to 
calculate average current (sub-section B) and the dc noise 
power (sub-section C) as functions of gate and source- 
drain voltage for two limiting cases: phonons equilibrated 
to the external world independent of the electron occupa- 
tion (phonon equilibration fast compared to dwell time 
of electron on molecule), and phonons uncoupled from 
the external world and responding only to on-dot elec- 
trons (phonon equilibration slow compared to dwell time 
of electron on dot). We find that in the case of phonons 
uncoupled from a bath, under certain bias conditions the 
phonon distribution can deviate far from equilibrium. 

Section-IV deals with the low temperature quantal 
regime of the phonon coupled dot where we use the 
Keldysh Green's function technique to calculate the dc- 
current and phonon distribution function to leading or- 
der in the ratio of (electron-phonon coupling) / (tunneling 
rate to the leads) for the two cases of slow and fast 
phonon equilibration rate. Section IV is split into sub- 
section A that re-introduces the problem, sub-section B 
that derives expressions for the exact eigenstates in the 
absence of electron-phonon coupling. These eigenstates 
form a convenient basis for carrying out the perturbative 
Keldysh calculation which is outlined in sub-section C. 
Sub-sections D and E present results for I-V for the two 
cases of slow and fast phonon equilibration rates. 

Finally Section V studies the crossover from low-T to 
high-T regimes and section VI is a summary of our work 
and its relation to already existing literature. 



II. MODEL, PARAMETERS AND FORMALISM 



Model 



with 



H D = en d + ^-n<j [n d - 1) + Xujo (b f + b) n d 

n 2 1 
+ -^ + -Kz 2 +uj b 1i b 
2M 2 



leads — ^ ^ ffc^Q, j^^CK,fe 
k,a=L,R 

Ht = 51 ti,k, a {z}a) ak d. 

a—L,R,k,(T,i—l,dg 



(2) 
(3) 

(4) 



where a labels the spin index. Fig.Q]shows the schematic 
of the energies considered. 

Here the number of electrons on the molecule, n d , is 
given by 



n d 



-l,d Q ,a 



(5) 



and the parameter U is the charging energy of the 
molecule. We have defined the zero-phonon state of the 
vibrational mode of the molecule to be the ground state 
when n d = 0, and we neglected anharmonicity in the lat- 
tice part of the Hamiltonian (such anharmonicity is of 
course induced by the electron-phonon coupling and an 
intrinsic anharmonicity could easily be added). 

The dot-lead coupling ti t k, a means that [H,n d ] ^ 0; in 
the absence of electron-phonon and many-body physics 
this implies an inverse lifetime for decay of an electron in 
state i on the dot into a state of energy e in lead a, 



(6) 



The density matrix p of the full Hamiltonian H obeys 
the equation of motion 



dp(t) 
dt 



= -i[H,p(t)} 



(7) 



The current / and the noise S through the lead a are 
given by 

(I a (t)) = Trp{t% (8) 

S a (t) = Trp(t)(SI(t)5I(0)+SI(0)5I(t)) (9) 



We consider the case of a molecule with a single level 
of degeneracy d g coupled to two leads, which we label as 
'left (L) ' and 'right (R)'. We suppose that the electrons 
are coupled to two different kinds of phonons; an internal 
vibrational mode of the molecule of frequency loq, which 
couples to the local charge, and a phonon mode labeled 
by a displacement operator z, that accounts for the os- 
cillations of the dot in an external confining potential 
of parabolicity K. This phonon mode does not couple 
directly to the charge on the dot, but results in an ex- 
plicit z dependence in the left and right tunneling matrix 
amplitudes ti^, a {z}- The full Hamiltonian is therefore, 



H = Ht 



(1) 



FIG. 1: Energy level diagram. pl,r represent the chemical 
potential in the left and right leads respectively, while e rep- 
resents the ground state of the singly occupied dot and the 
dashed lines indicate phonon excitations of the singly occu- 
pied dot. The solid line indicates the energy e" + U of the 
doubly occupied dot. 
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where SI(t) = I it) — (I) and the current operator through 
the a lead is given by 



(10) 



^ (*i,a(z)a*a,fc,<^ - h.c) 



N a = a a \ cr ak ,<? being the number operator for the 
a lead. 



B. Parameters and regimes 

The behaviour of the models we consider is specified 
by two important dimensionless parameters: the ratio of 
the temperature T to a typical decay time T, and the 
product of the dimensionaless coupling A and the ratio 
of the phonon frequency uo to T. For large values of 
T/r a classical rate equation analysis may be employed 
for all values of Xujq/T; this is the subject of Section III. 
For small values of T/T a quantal treatment is required. 
Section IV reports results for low T, obtained using per- 
turbative calculations valid for Xuiq/T < 1, while Sec- 
tion V treats the quantal-classical crossover, also in the 
Xloq/T < 1 regime. The quantal strong coupling regime 
(T/r < l,Aa>o/r > 1) is a challenging problem left for 
future reserach. The different regimes and the sections 
treating them are shown in Fig. [21 



C. Formalism 

The essential assumption in the theory of molecular 
devices is that the leads are in equilibrium independent 
of the state of the molecule. In order to implement this 
assumption it is often convenient to define a projected 
density matrix p s 



Ps = Tri eads {pit)} <g> pi ea( ls 



(11) 



T/r 



Classical Regime 
(Section III) 



Quantal- Classical 
Crossover (Section V) 



Perturbative Quantal 
Regime/Section IV) 



Strong Coupling 
Otiatihd Regime 
(Not studied here) 



FIG. 2: Overview of different regimes studied. 



where pleads is the density matrix of the left and right 
leads that are in thermal equilibrium at some specified 
chemical potentials (fiL and ijlr). The density matrix of 
the dot (also known as the reduced density matrix) is 
therefore given by 



Pd 



Tvi ea dsPs 



(12) 



and for our model corresponds to projecting the full den- 
sity matrix to the smaller subspace of the degrees of free- 
dom of the dot electrons and the two types of phonon 
modes. 

A complementary density matrix p t may also be de- 
fined such that 



Pt = P - Ps 
from which it follows that 



Tr. 



leads 



Pt = 



(13) 



(14) 



It is desirable to obtain the reduced density matrix of the 
dot pd because its diagonal components relate directly to 
the occupation probabilities of the various states of the 
dot, and various expectation values (such as the current 
and noise) may be expressed in terms of components of 

PD- 

In what follows we outline a general scheme for ob- 
taining the reduced density matrix 17 . The equation of 
motion for the full density matrix, Eq. [7] implies 

dp s i(t) (rp dpi(t) 

= I J ri eads 



>, i - ■ m . i Pleads 

at \ at 

= -i iTrieads[H t I it), pit(t)])® Pleads (15) 

dptijt) = dj Pl jt)~ p sI jt)) 
dt dt 

= -i[H tI (t),p sI (t)+p tI it)] 

+ iiTrieads[H t I,Ptl(t)}) <£> Pleads (16) 

We have used the interaction representation defined by 



Oj(t) 



- MH D +Hu 



0(t)e 



-i(H D +H u 



,)t 



In specifying the solution of Eqns. ^| and we 
require an initial condition. There are two common 
choices: (i) at the initial time, p corresponds to an equi- 
librium ensemble for H with /i£ = pn, (ii) at the initial 
time (ti) the dot and leads are decoupled (7J t = 0) so that 
p{t%) = Pd ® Pleads with p D and pi ea ds the equilibrium 
density matrices corresponding to the uncoupled prob- 
lem. We shall be interested in steady state, so we will 
take the initial time U — — oo. Further we shall be inter- 
ested in cases in which the leads are not strongly affected 
by the presence of the dot. In this case the boundary 
condition (ii) is most convenient, but we note that if or- 
thogonality effects or Luttinger liquid renormalizations 
of tunneling amplitudes are important then this choice 
may be less convenient. 

Choosing boundary condition (ii), Equation 1161 being 
a linear equation in pti(t) may be formally solved as 

/•OO 

pa(t) = -i dt / K I it,t')[H t iit , ),p sI it / )} (17) 
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where Ki(t,t') obeys the following operator equation, 



dKi(t,t')m 

dt 



+ i[H tI ,K I (t,t')* 



-i(Tri eacla [H t i,Ki(t,lf)»])®pi ead3 = S(t-lf) (18) 

where the symbol • denotes the operators acted on by 
Kj. The combination of the second and third terms on 
the l.h.s. of Eq. El correspond to processes in which the 
particle distribution in the leads differ from the equilib- 
rium one. When substituted back in Eq. El these terms 
produce correlations between the dot and lead variables. 
Processes corresponding to these correlations have been 
discussed in detail by Scholler and co-workers^, however 
their importance for the present problem is unclear. 
Substituting Eq. El m to Eq. El leads to 



dpDi{t) 
dt 



Tr leads / tf[H t i{t),K 1 (t,l!)[Ht J (l!) 1 p.itf) 

J — CO 

(19) 

Eq.EJis a generalized Master equation that under steady 
state conditions (setting the l.h.s. of Eq. EJto zero) can 
be written in the following form, 



(20) 



where (i\po\i) = P 1 , the probability of being in the i- 
th dot state. (Under steady state conditions off-diagonal 
elements of the density matrix vanishes except in the case 
of accidental degeneracies). Conservation of probability 
requires J^i Rj^i = 0; from which it follows that 







(21) 



The above expression for p t when substituted in the 
first line of Eq. El leads to the widely accepted Meir- 
WingreeniS expression for the current, 

II = J^dt' f deN(e)Im{e*< t - t ">G<(t,t>) + f L (e\25) 
e i<t-t')G R (t,t')} 

where N(e) is the density of states, and the Green's func- 
tions are defined as 

G R (t,f) = -i6{t-md{t)k L {m}J{^ k LW)m) 

G<(t,t') = i(dHt% L {z(t')}d(t)t k L{m}) (27) 

In general the second and third terms in Eq. 1181 im- 
ply Kj(t) — 9(t)k(t) where k ^ 1, suggesting the 
Meir-Wingreen formula is incomplete. However we sus- 
pect that in the present problem in which there is no 
backscattering and no "excitonic" interaction between 
dot and lead variables, correlations corresponding to 
i[H t i,.} - iTri eads [H t i, .} <8> pieads ^ (Eq. EJ) are ir- 
relevant: any such configuration simply propagates away 
from the lead and does not return. Some supporting 
evidence for this argument is provided by the exact diag- 
onalization results of Section III; however the issue war- 
rants further investigation. 

For the simpler case of left and right phonon indepen- 
dent tunneling amplitudes the Green's functions in Eg. 1251 
are simply the dot d-electron correlator. Wingreen et 
al4& showed that in this scenario and for the case of left 
and right tunneling amplitudes that are proportional to 
each other (but with non-trivial on-dot interactions), the 
current through a single resonant level simplifies to 



so that the quantities Rj—>i may be interpreted as the 
various in-scattering and out-scattering rates. 

A formal expression for the current may also be derived 
by using Eq. [HJand observing that Trlp s = 0, so that 

(/) = Trlp t (22) 

p CO 



Eons. PHI and 1221 are our basic results. Further analysis 
depends on the specifics of the system and of the approx- 
imation chosen. At temperature T, only times t < 1/T 
are relevent, so at sufficiently high T, we may approxi- 
mate Ki by its short time behavior, 



/ = 



de r L (e)T R (e) 
2ir r L (e)+T R (e 



(28) 



where A(e) = i(G R (e) - Gj(e)), with G R (t,t') = -i6{t- 
t')({d(t),dHt')}+). andr i/fl (e) - 2tt E k tl L/R S(e k -e). 

When the temperature is the highest energy scale in 
the problem, approximation of Eq. 1231 becomes exact, 
and further the Green's functions in Eq. El are replaced 
by their short time behaviour—. This shall be explic- 
itly demonstrated in the next section where we will de- 
rive rate equations for the on-dot probability distribution 
within the sequential tunneling regime. Following this, 
in Section IV we shall carry out the Meir-Wingreen pre- 
scription for calculating the current for general tempera- 
tures. 



Ki{t,t') -> 6(t-t') 



(23) 



Indeed making this substitution in Eq. El leads to the 
following expression for p t in the Heisenberg representa- 
tion, 



Pi 



-i /' dt' e-^n+^^-^Hup^t')] (24) 



i(H D +H Uads )(t-t') 



III. HIGH-TEMPERATURE APPROXIMATION: 
RATE EQUATIONS 

A. Formalism 

In this section we shall carry out a high-T analysis for 
the simple case where the only phonon mode the electron 
degrees of freedom couple to is the on-dot vibrational 
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mode. (The case of the phonons coupled to leads has 
been discussed is some detail by McCarthy et aliii within 
the high temperature approximation). 

It is often convenient to choose a representation which 
is diagonal in the dot degrees of freedom. In the present 
model this is achieved via a standard^ canonical trans- 
formation. Defining S = A (j>2i a d\ a di^ (b^ — b) and 

transforming all operators O via e Oe~ leads to a 
transformed Hamiltonian H' = H' dot + H' t + Hi eads with 



H' 



e'n d + Lo ttb + %n d (n d -l) (29) 

H.c) (30) 



H-t ^J2a=L,R,i^,aJ2p,<7\-^- a apa-di,a 



where the transformed phonon operator b = b — 
A J2i a *H a d i_<j, so that the phonon ground state depends 
on the dot occupancy. Moreover e' — e — X 2 uiq is the 'po- 
laron shift' in the energy for adding one electron to the 
molecule and the interaction parameter U is also renor- 
malized, but as we shall focus here onU — > oo we do not 
write the renormalization explicitly. The crucial phonon 
renormalization of the electron-lead coupling is given by 



X = cxp -A (&t - 6 



(31) 



The high-T approximation proceeds from Eq. 1241 by 
making the Markov approximation which involves replac- 
ing p s (t') in the integrand in Eq. HUby p s (t). After sub- 
stituting for p t in Eq. 1151 it is also convenient to replace 
J_ = i j which amounts to absorbing any level 
shifts induced by the dot-lead coupling into the bare val- 
ues of the parameters. Following this we obtain 



= -i[H D + H leads , Ps ] - | dt> (; 

[ Ht:e -i(H D +H leads )(t-t') [H t , Ps (t)}e i ( H ° +H ^X t - t "> 



On the assumption that orthogonality effects may be 
neglected, we may formally take the trace over the lead 
degrees of freedom and in the process arrive at coupled 
equations of motion for the various occupation proba- 
bilities of the dot. We outline the calculation in detail 
for one of the 4 terms that one gets on opening up the 
commutator in Eq. 1321 

Tr leads {dp ^ d fP°) (33) 
= -\^ 00 dt'Tr leads H t e^ H °+ H ^ t - t ') 

HtPleads ® p D e^+ H '^^ t - t "> +... 

Using the following relations 

TrieadsPleads = 1 (34) 
r leads {^PleadsO i a ki ap t k^j = S a ^Sk 1 ,fc 2 /(efc — p a ) 



Trt, 



Tr, 



leads {^Pleadsa a ,ki a \,k^j ^ ^a^^kxM (1 — f{ e k — Pa)) 



we obtain 



Xe- iH ^-^d] a X^p D e m ° (*-*') 

+ 4(l-/fe,a)e- iefc '" (i - t ' 3 
4xt e - iffD C*- t ')rf, CT Xp D e iffD ( i - t ') + ■ • • 



(35) 



We may now identify the probability of the dot being 
in a state with n electrons and m phonons as, 



P" 



(n, m\pjj\n, m) 



(36) 



Note that while po is always diagonal in the dot electron 
number, it may be off-diagonal in the phonon number 
(due to the presence of the X operators in Eq. I35|) . How- 
ever such terms are negligibly small in comparision to the 
components of po that are diagonal in phonon number, 
and are therefore neglected in our analysis. That is also 
the reason why we have dropped the first term in Eq. 1321 
in the next set of equations. 

We are now in a position to write rate (Master) equa- 
tions for the electron- phonon joint probabilities, which 
take the form 

^ = E 0)9 ' /« ((<? ~ <?') + U(n 1)) r« g ,P, ( r 1} (37) 
+ (1- fa(W -q)^ + Un))Tl q ,P^ +1) 
-(I- fa {{q - q') oJo + U(n - 1))) I^P™ 
-f*{tf-q)uo + Un) T a qlq P^ 

Note that in our notation the upper index in P™ always 
refers to the electron number and the lower index the 
phonon number, while f a (x) is short form for the Fermi 
function evaluated at f(x+e'—p a ), p a being the chemical 
potential of lead a. 

Thus within the high-T approximation the rate for 
going from an n electron and q phonon state on the 
dot to an n — 1 electron q' phonon state is R q Z^^~ = 
J2 a =L,R fa((q~ q') + U(n - 1)) T qql , where rj >9 rep- 
resents the transition rate involving hopping an electron 
from the dot to lead a and changing the phonon occu- 
pancy from q (measured relative to the ground state of 
H' D with occupancy n) to q' (measured relative to the 
ground state of H' D with occupancy n — 1) and is equal 
to the transition rate involving hopping an electron from 
the lead a to the dot and changing the phonon occupancy 
from q (measured relative to the ground state of H' D with 
occupancy n — 1) to q' (measured relative to the ground 
state of H' D with occupancy n). More explicitly 20 



r»,, q = r a \< q '\x\ q >\ 



(38) 



The matrix element can be computed by standard 
methods^; its absolute value | < q\X\q' > | 2 = Xj, n , 



qq' 
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is symmetric under interchange of q and q' and is 



X. 



q<q' 



k=0,q 



(-A 2 ) fe ( g !g'!) 1/2 Al«-«'le- x2 /2 



(k)\((q-k)\(k- 



(39) 



As interesting special cases, we write several lowest 
operators: 



Xqt, 



X U = (l - A 2 ) 
X 2 i =>/2A(l-£ 
A 22 - f 1 - 2A 2 + 4 



= -A 2 /2 



e -A»/2 



(40) 
(41) 
(42) 

(43) 



Observe that for certain values of A some of the matrix 
elements vanish. This unusual behavior is an interference 
phenomenon, which is slightly obscured by the notation. 
A state which has q phonons excited above the ground 
state of the system with n = electrons is a superposition 
(with varying sign) of many multiphonon states, when 
viewed in the basis which diagonalizes the n — 1 electron 
problem, and therefore the transition described by X qq i 
is really a superposition of many different transitions, 
which for some values of A may destructively interfere. In 
several recent papers^ the phonon renormalization of the 
molecule-lead coupling is apparently omitted, or treated 
in an average manner which neglects the q, q' dependent 
structure. 



B. I-V characteristics 

In this subsection we shall discuss the I-V characteris- 
tics obtained from the solution of the high temperature 
rate equations for two extreme cases. One is a scenario 
where the phonons are not coupled to a bath and their 
number changes only when electrons hop on and off the 
dot. The second case is when the phonons are strongly 
coupled to a bath, and are always forced to be in equi- 
librium. 

From Eq. we obtain 



(I) = Trp t (t)I L 



with 



=L/R 



it L 



£(< 



dX - d)xh 



(44) 



(45) 



Using Eq. [31 for pt and following the same procedure 
of tracing out the metal degrees of freedom, we arrive at 
the following expression for the current through the lead 
a in terms of the joint probability distribution functions, 

la = (2rf s - n)P»f a ((<?' -q)u; + Un) 1^(46) 

-(n + l)P q l+1 (1 - /„ {{q - q') u; + Un)) T«, q 



where the sum on n is from to (2d g — 1), 2d g being the 
maximum occupation of the dot. 

As an aside consider the simple case of spinless electron 
with no coupling to phonons with a non-degenerate dot 
level of energy e^. In this case the rate equations give us 
the following probability for a singly occupied level, 

pi = r L /(e rf - fi L ) + T R f(e d - ix R ) ^ 
El + E# 

while the current from Ea. 1461 is simply 

II = r7Ti% {/(ed - - /(ed - fiR)} (48) 

On comparing the above expression with the exact solu- 
tion for the current obtained by Wingreen et al. (Eq. I28|) 
we can explicitly see that the high-T approximation cor- 
responds to assuming the spectral function is a delta- 
function. 

We shall now discuss the opposite limit, of phonons 
equilibrated to an independent heat bath, assumed to 
be at the same temperature as the leads. To implement 
this we force the probability distributions on the right 
hand side of Eq.|37|to have the phonon-equilibrium form 
P r q l = P^e- quj °f T {l - e-^l T ). In the U -> 00 limit this 
ansatz implies that the probability P° that the dot is 
empty is given by, 



P° = 



E-na p -qu a /T s 
a.q.q' q,q' J a,q,q' 



Ea.q.q* ^lq>e-^o/T faqql+T a^ e - qUo /T fa ^ 

_(49) 

where f a q q , = 1 - f a ((q - q') u ), fa,q,q> = 1 - /„,,,,' 
while, P 1 = 1 - P°. 

In general for both equilibrated and unequilibrated 
cases the rate equations may be written in the matrix 
form 



P = MP 



(50) 



Therefore under steady state conditions (P n = 0), the 
problem reduces to finding the eigenvector correspond- 
ing to the zero eigenvalue of the matrix M. We do this 
numerically. From these solutions we have computed 
the current. Representative results are shown in Fig. |3 
which plots the low-T current as a function of V s d for 
two gate voltages: V g = (pl = —j-iR, upper panel) and 



V a 



5* (Mfl = 0j lower panel), for both equilibrated 



and unequilibrated phonons. (Note that the calculation 
is for large values of U which correspond to negligible 
double occupation probabilty for the dot electrons). 

Steps (broadened by T) in the current associated 
with "phonon side-bands" are observed when the 
source-drain voltage passes through an integer multiple 
of the phonon frequency However, in the opposite 'linear 
response' limit V s d — > (not shown), as V g is varied we 
find just one main step in the I — U-curve, as V„ passes 
through 0, and only very tiny structures (vanishing as 
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FIG. 3: Current (J) vs source-drain voltage V s a for coupling 
constant A = 1.0 ujq = 1 and T = 0.05. Upper panel is for 
V g — 0.0, while lower panel is for V g = V s d/2, lir = 0. I is in 
units of ^f- 

h 



FIG. 4: Ratio of differential conductance (G) peak heights 
for equilibrated and unequilibrated phonons and two different 
coupling strengths and fiL = — fiR- The points for A = 0.5 
(open symbols) have been multiplied by 10. 



e _w °/ T , which is the probability of the dot being empty 
with one phonon excited) when V g is a non-zero multiple 
of the phonon frequency. This result appears to differ 
from that stated by other authors who find phonon 
side bands as V g is swept at V s d — > 0. The authors of All 
apparently neglected the fact that 

the phonon side-bands "float" i.e., shift with the fermi 
level as V g is changed. 

Fig. reveals on first sight an apparently surprising 
result: for symmetric bias (V g = 0) and for the coupling 
considered, the current is larger for equilibrated phonons 
than for the unequilibrated case, whereas for the strongly 
asymmetric case (lir = 0), the opposite is true. This is 
surprising because one expects that in the unequilibrated 
case the phonons arrange themselves so as to maximize 
the current. To gain more insight into this phenomenon 
we have calculated the dependence of the ratio of currents 
for unequilibrated and equilibrated phonons on the cou- 
pling A for different degrees of bias asymmetry. We find 
that except for lir — 0.0 (the most asymmetric case) a 
minimum in the ratio occurs for a A~ 1. This behaviour 
may be traced back to Eg. 14 11451 which reveal that higher 
order "diagonal" (n phonon- n phonon) matrix elements 
vanish for a A~ 1. 

The steps in current may be conveniently parameter- 
ized by the height (or the area, as the width is sim- 
ply proportional to T) of the corresponding peaks G max 
in the differential conductance G = dl/dV. Ratios of 
peak heights (or areas) provide a convenient experimental 
measure of whether the phonons are in equilibrium. At 
low T, the equilibrium phonon distribution corresponds 
to occupancy only of the n — phonon state, so the n-th 
side band involves a transition from the phonon to the 
n phonon state. Therefore the ratios of the peak heights 



-St 



2 4 6 



2 4 6 8 10 



FIG. 5: Phonon probability distributions for two differ- 
ent electron-phonon coupling constants calculated for fit = 

—fXR = 2u)Q. 



Eans. |4*U1 ETUI imply that if lil — —lir and T«w , 



Gn 



G° 



\x 



/),() 



equil 



2\X, 



oo 



2(n\) 



(51) 



or areas are controlled by ratios of \X„ 



In particular 



Note that if ill, [ir, then [i dependent changes in the 
occupation probabilities lead to additional, and not sim- 
ply characterized n dependence. 

Deviations from this pattern imply non-equilibrium 
phonons. As illustration we display in Fig. 0] G max val- 
ues (normalized to the zero frequency peak) for equi- 
librium and non-equilibrium phonons and a weak and 
strong electron phonon coupling. One sees that in the 
non-equilibrium case the peak heights display a non- 
systematic dependence on electron-phonon coupling and 
peak index, but that in general measurements of the 
n = 1 and n = 2 peaks reveal the effect clearly. 

It is also of interest to consider how far out of equi- 
librium the phonon distribution may be driven. Fig. [S] 
shows the phonon occupation probabilities for weak and 
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n n V 



FIG. 6: Phonon probability distributions for weak electron- 
phonon coupling (A = 0.5) calculated for /il = 2u>o,p,R ~ 0-0 
(left panel) and pl = 5o»q, \ir = 0.0 (right panel). Note 
the saturation in the probability distribution function which 
is also much closer to equilibrium. This situation is quite 
different for the same coupling constant and symmetric bias 
(see Fig. 03 • 



strong electron-phonon coupling and V g = 0. One sees 
immediately that the phonon distribution function is far- 
ther from equilibrium for weak couplings than for strong 
couplings. We associate this effect to the strong A de- 
pendence of operators Xq u , (Eq. I40|) which allows the 
system at large A to "jump down" from a highly excited 
state to one of low phonon occupancy. The deviation 
from equilibrium is largest for V g = for similar reasons. 
Fig. illustrates the scenario of non-zero gate voltage 
or asymmetric bias conditions pl = V s( i,Pr = 0.0. Here 
the phonon distribution for weak coupling saturates with 
bias to a value which is closer to its equilibrium distri- 
bution. As we shall show in Section IV, this gate voltage 
dependence of the non-equilibrium phonon distribution 
function is recovered in the quantal regime r > Tas well. 
Fig. [7| is the average phonon number N p h — J2 n m n ^n 
for moderate electron-phonon coupling. The steps in 
Nph — V s d observed here coincide with the steps in I- 
V and correspond to sequential (direct) tunneling. This 
is to be contrasted with the quantal regime (Section IV, 
Fig. 115(1 where N p h increases continously with bias due 
to higher order cotunneling processes. 



C. DC Noise characteristics 

Another important spectroscopic tool that is sensitive 
to the details of the electron-phonon coupling and to the 
phonons distribution is the current noise. In this sub- 
section we outline the calculation of the dc current noise 
within the high-T approximation. Quite generally, the 
current noise S ll (t) through the left lead is given by the 



FIG. 7: Average phonon number under symmetric bias con- 
ditions, hl = ~fJ,R , electron-phonon coupling A = 1.0 and 
uj = 20T. 



following correlation functioniiSl 

S LL (t) = ^Trp{I L (t)I L (0) + I L (0)I L (t)} - (TrpI L f 

(52) 

Using the fact that I L (t) = e lHt I L (0)e~ tHt , the above 
expression for noise may be rewritten as 

S LL (t) = ~Tr\(t)I L (0) (53) 

where 

A(i) = e- mt {(I L - (I L )) p + p(I L - (I L ))}e lHt (54) 

Since we are calculating the correlation of the same phys- 
ical quantity at two different times, we expect S(t) = 
S(—t). Therefore we shall explicitly calculate S(t) for 
positive times, for which we need to calculate the causal 
function X(t) that obeys the equation of motion 

^ = -t[H,\(t)]+S(t){6I LP + p5I L },t>0(55) 
= 0, t < 

where 81 l = Il - (II)- Note that Tr A(0) = 0, and we 
expect X(t) to be traceless at all times. Also from charge 
conservation it follows that the noise across the left and 
right leads are equal = Sim(t)). 

We now solve the equation of motion for A by decom- 
posing A = A s + At, where A s = Xd <8> Xieads is diagonal in 
the dot and lead variables, while A* is off-diagonal. After 
doing a similar decomposition for the density matrices 
p = p s + Pt, the equation of motion for A (Eq. I55|l may 
be re-written in a manner similar to that done for the 



9 



equations of motion for p in the previous section 

-%[H t , X t ] + 5{t) (I Pt + p t I - 2(I) Ps ) (56) 



dt 
dt 



— —i[HD + Hl ea ds, At] — i[H t , A s ] 
+ S(t)(I Ps +p s I -2(I) Pt ) 



(57) 



Note that we have dropped the [Hu+H m , X s ] term for the 
same reasons as before, namely that the matrix elements 
of this term is off-diagonal in the phonon number and is 
very small in comparision to the components of A s that 
are diagonal in phonon number. 

The solution for At from Eg. 1571 is given by, 

\ t {t)=-if dt'e-^+^^-^lH^Xsit')} (58) 

J — CO 

e i(H D +H leads )(t-t') + Q, t \ e -i(H D +H lcads )t 

(J(0)p s (0) + p s (0)7(0) - 2(l)pt(0)) e^ +ff —>' 

Substituting Eq.|5Hlm Eq.lSrJland going through the usual 
steps of extending the upper range of the integral to in- 
finity, and making the Markov approximation (which in- 
volves pulling A out of the time integral) , we arrive at the 
following matrix expression for the equation of motion for 
A, 



dX, 



dt 



= MX D +S(t)h, t > 
= 0, t < 



(59) 



The matrix M is the same that enters in the equation 
of motion for pp. In arriving at the expression for the 
vector h, we have approximated oscillating factors such 



iet 

structure 



as e"" ~ 2it8{t)8{e). Following this the vector h has the 



(n + i)p; +1 4 



r>L,n,n+l 



9 -l-q' 
X/,n+l,n 



(60) 



where as before the sum on n is from to 2d g — 1, d g 
being the number of degenerate levels not counting spin. 
Note that we are using the following short-hand notation, 



R 



a,n,n+l 
q,q> 



K Q',g 



f a ((q'-q)cjo + Un)ri q , (61) 
= (1- fa ((q'-q) wo + Un))Tl q , (62) 



Moreover by using Eq.0|j]it is easy to check that h (whose 
components are given in Eg. I60|) is traceless. 

Now we shall rewrite the Eq. 1531 explicitly in terms of 
the components of Xd and the steady state probabilities 
P® and P„ (components of reduced density matrix pd)- 
After some algebra one finds that 

Sll (t) = ^Tr m , D I L X t , t > 

= -{Il? + \{S lL {t) + S 2L {t)) (63) 




FIG. 8: DC Noise for A = 0.2 and under symmetric bias 
conditions pl = — PR = V s d/2. V s d is in units of cjq = 20T. 
Sdc is in units of 



where S L L(-t) = S L h{t), and 

SiUt) =E q ., ql .M-n)X-(t)R^' n+1 
-(n+l)A«+H<)<'; +1 '" 

and 



S 2L (t) = 25(t) E,,^,„(2d fl - n)P?(t)R 
+ (n+l)P«+Ht)<J +1 '" 



(64) 



(65) 



Note our notation A™ = (n,q\X D \n,q) etc. where A D = 
TrieadsXs- Also note that the S(t) in the expression 
for S2L(t) again arises from replacing oscillating factors 
e let ~ 2ir5(t)S(e). 

We find it convenient to perform the following shift 
of variables X D -> X D + 2(1 L ) P. This shift of variables 
does not affect the equations of motion Eq. EH1 because 
MP = in steady state. However this shift of variables 
cancels the (II) 2 term in Eci.l6"31 

Collecting all the terms we arrive at the following ex- 
pression for the DC noise power 



Sdc — 2 



dtS LL (t) 



(66) 



S, 



= 2£^,„(2d 9 -n)(P 9 , 

-(n+l)(-P^ + ~X^(0))KT 1 ' n 



where the components of X(lu = 0) are obtained from 
solving the matrix equation 



P(w = 0) 



(68) 
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Note that the matrix M has the property that ^ My = 
0, so that one of its eigenvalues is 0. However h being 
traceless, M _1 /i is well defined. 

Let us now look at the simple case of no electron- 
phonon coupling and spinless electrons. In that case 
the matrix M acqires the simple 2x2 form (note that 
/ a = l/(l + e^«)/ T ) 



, f _ - {TlU + T R f R ) T L (1 - f L ) + T R (1 - f R ) 
1 T L f L +r R f R -(T L (l-f L )+T R (l~f R )) 



while the components of A(w = 0) are given by 

h 1 



A^-A^- 



Tl + Tr 



(69) 



(70) 




where 



h 1 = 2T L P°f L - 2(I L )P 1 



(71) 



The full expression for the dc-noise is 



S dc = 2T L (P°f L + P 1 (l-f L ))-4 



;P°h (72) 



T t +r, 



MP 1 



We may now derive expressions for the dc noise in two 
limits. The first one is in the linear response regime [II = 
HR so that P° = 1 — /l and (II) = 0. In that case, 



FIG. 9: DC Noise for 
conditions /il = — = 
Sdc is in units of ^-r-. 



A = 1.0 and under symmetric bias 
Vsd/2. V s d is in units of u>o = 20T. 



IV. QUANTUM THEORY OF TRANSPORT 
THROUGH A PHONON-COUPLED DOT 

A. Overview 



i-L+i-R 



(73) 



which is the result expected from the fluctuation dissipa- 
tion theorem Sdc = 4TG. 

The other limit is = —fiR = eV/2 > T so that by 
using the expression P 1 = ri -^^p^ a and setting to zero 
combinations such as /l(1 — /l) we obtain 



S{w = 0) 



F 2 +F 2 

i L ' R 

'(Xl + Tr)- 



(74) 



where (I L ) = ^ l l ^ r . The above expression gives the 
standard shot noise result when <C Fr — 

The results for the noise for the case of phonons in equi- 
librium and the opposite case of phonons not coupled to 
any heat bath are illustrated in figures |H] (weak coupling) 
and[5] (strong coupling). The difference between the equi- 
librated and unequilibrated cases is more dramatic for 
smaller A. For A = 0.2 while only one phonon side band 
is seen for the equilibrated case, very sharp phonon side 
bands are seen for the entire bias range when the phonon 
distribution is far out of equilibrium, with certain side 
bands appearing as peaks (rather than steps) associated 
with the suppression of noise. 



In this section we present a fully quantum mechanical 
treatment of the simple limit of the model considered in 
previous sections. In order to carry out the calculation 
we restrict attention to a single non-degenerate level with 
no onsite coulomb U, and to low orders in perturbation 
theory in the clcctron-phonon coupling. The results shed 
light on the relation between the Green's function for- 
malism natural in the quantal treatment, and the density 
matrix formalism natural in classical problems, and elu- 
cidate the quantal to classical crossover. For the reader's 
convenience we reproduce here the limit of Eans. l2l3l and 
0] which we study. 



Hei = e d^d+ £kal a ak,a 

k,a=LM 

Q>k,a ~r fl.CJ 



k.a—L^R 



H p h 
Hel-ph 



(75) 

(76) 
(77) 
(78) 



Note that we have neglected the spin index, which may 
be restored in the final expressions for the current by 
simply multiplying by a factor of 2. 

We use two methods to analyze the above Hamilto- 
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nian. One is the Keldysh Green's function method and 
the other is an explicit construction of the eigenstates 
and thus the density matrix. We note that the Green's 
function method directly computes the expectation value 
of operators at different times, bypassing the explicit con- 
struction of the density matrix. However equal time mul- 
tiparticle Green's functions correspond to moments of the 
density matrix, and permit in principle its reconstruc- 
tion. 

The rest of the section is divided as follows. In sub- 
section B we provide the exact solution in the absence of 
coupling to phonons. The eigenstates of this system will 
form the basis for the perturbative calculations that fol- 
low. In subsection C we outline the Keldysh calculation 
while in subsection D we present results for the special 
case when the phonon distribution is always in its ground 
state. In section E and F we generalize to the case when 
the phonons are allowed to deviate from equilibrium and 
also supplement our results obtained from the Keldysh 
Green's function technique by a perturbative calculation 
for the phonon density matrix that allows us to obtain 
the out of equilibrium phonon distribution function. 



B. Non-interacting dot and exact expressions for 
the current for an interacting dot 

Following standard methods 2 ^, the exact eigenstates 
of H e i fEa .175(1 can be easily obtained. The Hamiltonian 
after diagonalization is 



given by, 



E 

k,a—L,R 



e kal a a k 



(79) 



Note that while a^^L/R in Eg. 1751 refer to states that live 
only on the left /right lead, an exact eigenstate of H e i has 
non-vanishing amplitude in both leads. We write the ex- 
act solution in a scattering state basis in which afe !a =L/i? 
refers to a running wave incident from the left /right lead 
with a certain amplitude of getting reflected back to the 
starting lead and a corresponding amplitude for trans- 
mission to the other side. The label L/R now refers to 
the lead from which the particle is incident, and there- 
fore determines the distribution function describing oc- 
cupancy of the states. The dot and metal electron cre- 
ation/annihilation operators are related to these exact 
(scattering) eigenstates as follows 



0<k,a — ^ Vka,k'bClk',b 
k',b=L,R 

d = ^ Vk^k,a 
k,a— L,R. 



(80) 
(81) 



The coefficients rj and v in the above set of equations 
obtained from ensuring proper commutation relations are 



11ka,k'b = 5k,k'5a,b 



tkaVk'b 



Vka 



£ka ~ Cfe'6 + iS 
tka 



L k'h 



(82) 
(83) 



'b e ka -e h , b -i5 

We identify tunneling rates to the left and right leads by 



rz,/.R.(e) = 27r>Ji| L/fl i(e-efc 



k 



(84) 



The current through the left lead is given by 
dN L 



h 



dt 



t kL (d^akL - h.c.) (85) 



(h) = 2|^t fei /m(d 1 a kL ) 

k 

and charge conservation requires II — —Ir- Plugging 
in the expressions for d and akL in terms of the exact 
eigenstates aka into the expression for the current (Eq. 
1851 we obtain 



t L Im (ykLte[)Vku*( a kl, a a M,p) 



k : ki : k2',a,p — L,R 

(86) 

In the non-interacting limit the expectation value of the 
exact eigenstates is 



( a k ia a k2,p) = S a ,f3S kl ,k 2 f( e k - M/S) 



(87) 



where f(x) is the Fermi-Dirac distribution function. Sub- 
stituting this in Eq. 1861 we obtain 

/ t \ e f de T L (e)r fi (e) 

= h J ^ r £ ( e ) + r a ( e ) {/(g " M£) " /(e " Ma)}i4(e) 



where A(e) is the spectral density of the dot and is given 
by 



A(e) = 



( e - eo -£'W) 2 + (^) 2 



(89) 



Note that E'(e) is related to T(e) = T L (e) + T R (e) by 
the usual Kramers-Kronig relation. For simplicity in our 
subsequent computation we assume energy independent 
density of states and tunneling amplitude tL,R, so that 
£' = and Tl.h are constants. Ea.l88lagrees with the ex- 
pression for the current that was derived by Wingreen et 
ali& employing the Keldysh non-equilibrium technique. 

In the presence of electron-phonon interactions, the full 
Hamiltonian (now including Ea. 1771 and En .178(1 takes the 
following form in the basis of exact eigenstates of the non- 
interacting system, 



k,a=L,R 

+ \u; a (b + b r ) 



(90) 



E V k,a v k<t,u\ a Oik'b 
k,k' ,a,6— L,i? 
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A standard method for studying the nonequilibrium 
problem posed by H is to define retarded Green's func- 
tion for the electrons on the dot 

G R (t,t')=-i6(t-t')({d(t),S (<')}+) (91) 
= -i9(t - t>) J2 ka k , b v ka ^, b {{a ka (t),4, b (t')} + )(92) 

and the Keldysh Green's function for the dot electron, 

Gf(M') = <*t (*')}_) (93) 

^ a ^ 6 ({a fca (t),4 f ,(t')}-)(94) 



; E 

ka,k'b 



Similarly for the phonons we define the corresponding 
retarded and Keldysh Green's functions, 

D R (t,t') = -i9(t-t')({b(t) + tf(t),b(t') + tf(t')}(£ip) 
D K (t,t') = -i({fo(t) + fe+(t),fo(i') + fo t (i')}+> (96) 

Note that the phonon (electron) Keldysh propagators 
in the equal time limit are iD K (t,t) = 2(1 + 2(6 t 6)) 
(iG^(t,t) = 1 — 2(d}d)) and are directly related to the 
average phonon (electron) number and therefore corre- 
spond to the first moments of the density matrix. Higher 
order equal time correlators ((b(t)tf (t)) n ) give higher mo- 
ments of the density matrix, enabling in principle the full 
reduced density matrix pjj to be reconstructed. 

The retarded dot Green's function for a single resonant 
level in the absence of phonons can be easily obtained by 
using Eq. 1571 

G R {t,t') = g R (h,t 2 ) = Wa\ 2 G R aM {ti,h) (97) 

ka 

where G R aM (t,t/) = -iO(t - t')({a ka (t), 4a (*')}+> = 
—i9{t — t')e~ ltk ( t ~ t '. It is now easy to see that in 
Fourier space the retarded Green's function for the non- 
interacting dot has the familiar form 



E 



Wka\ 



ka "2 



(98) 



In a similar manner the Keldysh Green's function in the 
absence of phonons is found to be 



g K {uo) = 2TTiJ2Wka\ 2 5(LJ-e k ){2f(e k - n a )-l} (99) 
ka 

= T L (1 - 2}{u - /i L )) + Y R (1 - 2f{uj - m )) 

(UJ 6 ) 2 + £ 

where f(x) — l/(exp(^) + 1) denotes the Fermi distri- 
bution function. 

Moreover in the non- interacting limit, the retarded 
and Keldysh phonon Green's functions (D R ' K ) defined 
in Eg . 1951 and 1961 take the following form in Fourier space 



D R (co) 



2U! 



uj 2 — u)q + idsgn(uj) 

U} 



(100) 



In order to calculate the current for the case when 
electron-phonon interactions are present, we shall use 
the result derived by Wingreen et at, namely that 
the current through an interacting dot is still given by 
the expression Eq. |5SJ but with the spectral density 
A(e) — 2iIm{G R (e)}, where G R is the d-electron re- 
tarded Green's function calculated under appropriate 
nonequilibrium conditions and with respect to the full 
Hamiltonian H. We shall carry out this prescription for 
calculating the spectral density in the next section. 



C. Keldysh Greens function method: perturbative 
analysis 

In the presence of non-zero electron-phonon coupling, 
the Dyson's equations we wish to solve may be written 
in the following compact form in 2 x 2 Keldysh space 22 



D 



n 



(102) 
(103) 



where 



^ d - 1 Gi 



is the local dot Green's function, and 



S = 



o iy 



Dq («) = -2ni{5(Lu + lo ) + 5{lo - uj )} coth ^ (101) 



is the electron self-energy due to electron-phonon inter- 
actions. A similar matrix structure for the phonon prop- 
agator D and polarization II in terms of retarded, ad- 
vanced (D R / A ,H RA ) and Keldysh (D K ,IL K ) components 
also exists. The non-interacting Green's functions gd and 
Do have components that have been explicitly calculated 
in Eqn. EFJ ESI HUHI and fTTTTl Note that the temperature 
enters explicitly via the bare electron and phonon Green's 
function. 

We analyze the equations perturbatively. The expan- 
sion parameter is ^j? 1 , and the leading non-trivial S 
and II are represented by the diagrams in Fig. ^| We 
write the perturbative expansion in the usual self-energy 
language, but we note that in contrast to the conven- 
tional band-electron case crossed diagrams for the elec- 
tron Green's function are not small relative to uncrossed 
diagrams, because the Green's function lacks the pole 
structure found in the translation invariant case. Our 
results for the electron Green's function and electron ki- 
netic equation should be understood to be perturbative 
in A. 

To leading nontrivial order it is sufficient to calculate 
the phonon self-energy II using the bare electron Green's 
function gd, but a correct calculation of S requires the 
use of the full D. 

The retarded phonon self energy (Fig. I10|l becomes 

U R (t,t>) = - l 2^l {g ^t,t') 9 K {t\t)+g K {t,t')g A {t\t)} 

(104) 
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while the Keldysh phonon self-energy is 



n *(t,f) =zi^{ g R {t ,t/)9 A (t',t) (105) 

+g A (t,lf)g R (t',t) + g K (t,tf)g K (t',t)} 

Going into Fourier space and using Eans. and 1991 we 
obtain the following expression for the real and imaginary 
parts of tl R {uj) and fl K (Lo), 

u R (u) = -A 2 {r L Ti( AiL ,^) + r L T 1 ( ML ,- w )(io6) 

+ T R Ti(jj,R,u)+T R Ti(fi R ,-u})} 
n? m (cj) = i\ 2 {r L T 2 (n L ,w)-T L T 2 (ii L ,-u) (107) 
+ ^rT 2 (hr,uj) -T r T 2 {hr,-uj)} 



ElM = ^coth(^) {T 2 {^ l ,lu) ~ T 2 {u L ,~u;))(W8) 



+ 



2r 
rir, 



+-^coth(jfa) (T 2 (u R ,u) - T 2 (fj, R , -u)) 
K C0t h{^pSL) (T 2 (u l ,lu) - T 2 (ur,-uj)) 



IT 

m th( ) (T 2 (ur : , w) - T 2 (u L , - w) ) 

where we define the following integrals 

dw 2 (w + w 2 - e ){l - 2f(u 2 - u)} 



Ti (u, uj) 



and 



2^ ((wa - eo) 2 + ^)((w + cj 2 - e ) 2 + 

(109) 



1 - 2/(a, 2 - ,l) 



y 




FIG. 10: Diagrams that correspond to the leading contribu- 
tion to the electron self-energy (E) and the phonon polariza- 
tion (n). 



to obtain, 
D*(w) 



2uj 



T ( uj) = ~[ —L 

2(filUJ> 2 J 2tt (( W2 -e ) 2 + ^)((u; + a, 2 -eo) 2 + f) 

(110) 

Analytic expressions for Ti(fj,,u>) and T 2 (/i,uj) may be 
obtained at zero temperature, and are given in Appendix 
A. 

Note that the combination Ti(/Z£,,r£,,a;) + 
Ti(fJ,L,T L ,-uj) is symmetric, while T 2 ((i l ,T l ,lu) - 
T 2 (pl,T w) is asymmetric with respect to u>. As 
a result, for all combination of couplings and applied 
voltages 



uj 2 — lUq + iSsgn(uj) — 2u>qII r (lu) 



- ffu A -u 



(116) 



(117) 



A" 



n£(-w) = n?» 

n,m(-w) = -n£,M 
n K (-w) = n K (o;) 



(in) 

(112) 
(113) 



Note that our parameterization in terms of f% is such 

that at equilibrium f^ h (x) — coth(x) = 1 + 2ub(x), ur 
being the Bose distribution function. 

Besides T, R it is also useful to evaluate the Keldysh 
component of the self-energy Y. K that we will later use in 
order to calculate the distribution function of the dot. To 
the perturbative order to which we work, Yj K is related 
to the non-interacting Green's function of the dot and 
the phonon Green's function as follows 



The retarded electron self-energy represented by the 
diagram in Fig. 1101 is 

E*( t) if) = l -^{g R {t, t')D K {t, if) + g K {t, t')D R (t, t')} 

(H4) 

Once we know the components of the polarization matrix 
II, we may use the Dyson equations 11031 together with 
the following parameterization for D K — 



D K = D R f7 h - h Ph D A 



(115) 



l\ 2 0Jn 



= g R (t,t')D R (t,t') 



(118) 



+ g A (t,t')D A (t,t') + g K (t,lf)D K (t,if) 



Now that we know the retarded electron self-energy, 
we can use Eq.|SHlto calculate the current which is given 
by, 

I(^ l ,ur) = ^ (ftfa) J % - u L ) - f{u - ur)) 

(w-e -S« (w,ML,MH)) 2 + (5-Sf m (w,ML,MR))" 1 ' 
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FIG. 11: Zero bias conductance (in units of 5-) as a function 
of level energy (or gate voltage) with (ujo = 2F, solid line) 
and without coupling to a vibrational mode (dotted line). 




In the next two subsections we present our results for 
the two extreme cases of phonons strongly coupled to a 
heat bath, and therefore always in equilibrium, and of 
phonons uncoupled to an external environment. 



D. Results : Equilibrated phonons 

In this sub-section we specialize to the case where the 
phonons are always in their ground state so that the 
phonon Green's functions D R and D K in equations 11141 
and 11181 are always calculated under the condition that 
Ml = Mi?.; while all the non-equilibrium effects are in- 
cluded in the electron Green's functions. An important 
property of S^(ij) is that at equilibrium it is zero for 
u) = 0, and from this and Eg. I119l it immediately follows 
that the zero bias conductance even in the presence of 
electron-phonon coupling has the form 

W-V-TwStf (120> 

for the case where the two lead chemical potentials in 
addition to being equal to each other are also aligned 
with the dot level e^. Fig. ^] shows the gate voltage 
dependence of the zero bias conductance for symmetric 
broadening and two values of A. 

Fig. 1121 presents our result for the current (conduc- 
tance) for the equilibrated phonon case at zero tempera- 
ture. The top panel is the difference between the current 
with and without electron-phonon coupling for two differ- 
ent bias conditions. Note that for asymmetrically applied 
biases, the current with phonons can take a value larger 
than that in the absence of phonons. This is due to a shift 
in the center position of the spectral density (the total 
area under the spectral density being conserved). The 



FIG. 12: Zero temperature current (in units of 4p and mea- 
sured relative to no phonon current 7o [upper panel] ) and 
conductance [lower panel] (in units of ^-) under asymmetric 
(fiR = 0,flL = Vsd) and symmetric (fi L = -fiR = V ad /2) 
bias conditions and for equilibrated phonons. The phonon 
frequency u>o = 2.0r and Tl = Tr = 0.5r. Moreover the 
electron-phonon dimensionless coupling strength is A = 0.25. 
Note that under conditions of asymmetric bias, the current 
can saturate to a value larger than that for the device without 
phonons (A = 0). 



lower panel is conductance for the same bias conditions. 
The conductance for the symmetric bias (Mi = — Mi? ) 
case shows two features, one at V s d/u>o = 1-0 and the 
other is a much broader feature at V s d/u>o — 2-0. While 
the former corresponds to the onset of in-elastic cotun- 
neling, the latter corresponds to the onset of sequen- 
tial tunneling. (Under the asymmetric bias condition of 
Hl = Vsd j Mi? = one observes only sequential tunnel- 
ing). 

The transition from cotunncling dominant current to 
sequential tunneling dominant current can be understood 
by studying how the imaginary part of the electron self- 
energy due to interaction with phonons (S^J evolves 
with bias. Fig. IT31 shows T,f m {u) for the symmetric bias 
case hl — Vsd/2 = — Mi? = w o/2 and for two differ- 
ent values of T/loq. For simplicity we have considered 
the case where = T/j. Under symmetric bias , 
increases rapidly for \u>\ > \u>q — V s d/2\, while in calcu- 
lating the current the spectral density is integrated from 
— V s d/2 to V s d l '2 (see Ea. lll9|) . Therefore clearly there is 
a threshold at V s d = ujq when the £j m (w = <^o/2) jumps 
by cx ii^pa , and this corresponds to the onset of inelas- 
tic cotunneling. As the voltage is increased further, the 
range of integration also increases to finally include the 
Lorentizian broadening centered around ui = u>o, and this 
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FIG. 13: Imaginary part of the electron self-energy (in units 
of r) due to phonons for the symmetric bias condition /il = 
—fiR — Vsd/2 — uio/2, and Xluo = IT, level energy eo = 0. 
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FIG. 14: Cotunneling dl/dV (in units of for bias condi- 
tions ul = —fin. = Vsd/Z, level energy eo = 4.0r, u>o = 2r , 
Xloo = ir and for phonons not coupled to any heat bath. 



corresponds to the onset of sequential tunneling. 

As r is made smaller, the size of the step in Y> R m (u) = 
u>o/2) decreases, thus decreasing the cotunneling con- 
tribution to the conductance. Moreover the lorentzian 
broadening in the self-energy around ojq also becomes 
narrower, so that the sequential tunneling peak in the 
conductance also gets sharper. 

Fig. 1 141 shows the temperature dependence of inelastic 
cotunneling under conditions of no coupling to a heat 
bath (unequilibrated phonons). Calculations for this 
have been performed in the regime eo > Mi > PR so that 
the current is entirely due to cotunneling. (The current 
under resonant conditions for this unequilibrated phonon 
case is discussed in detail in the next section). It is clear 
from Fig.^Jthat inelastic cotunneling shows up as a step 
in dl/dV that gets rounded very rapidly with increasing 
temperature. 



E. Results: Unequilibrated phonons 

For out of equilibrium conditions, we may derive a 
quantum Boltzmann equation for the mean phonon num- 
ber, which for weak electron-phonon couplings is iden- 
tified as (N ph ) = - 1+ ^=^\ (This is from using 
Eq. 11151 and the fact that for weak couplings D R m is al- 
most a delta function at the phonon frequency). There- 
fore Eg . II 1 71 rewritten under steady state conditions and 



at an on-shell frequency has the form 

= (N ph )(U K (LJo) + U R (Lj )-^ A ^o)) (121) 
- (1 + (N ph )) {n K (oj ) - (IF>o) - ieVo))) 

From this the phonon outscattering rate may be iden- 
tified as {n K (uj ) + n R (uj ) - U a (uj )}, while the in- 
scattering rate may be identified as {IL K (ujq) — (H r (ujo) — 
H a (u)q))}. (Nph) is plotted in Fig.EDfor a variety of bias 
conditions. The results here are similar to what was ob- 
served in the high-T classical calculation, namely that 
for bias conditions under which the dot is half filled or 
close to it, the phonons tend to go far out of equilib- 
rium. When the phonons deviate considerably from their 
ground state, the corrections to the electron-self energy 
become comparable to F and one is no longer within the 
perturbative regime. Therefore the results for the cur- 
rent and conductance that we present here (Fig. I16fl are 
for the case of ul = V s d,PR = 0.0, when the phonons 
acquire a steady state distribution at large biases that 
is not far from its equilibrium value and a perturbative 
approximation is valid. 

The upper panel of Fig. ^j]plots the difference between 
the current with and without phonons for the asymmetric 
bias case, and for comparision this is plotted for both 
equilibrated and unequilibrated phonons, while the lower 
panel is the conductance peak corresponding to the first 
phonon side-band. Again within perturbation theory the 
differences between these two cases is not significant. 
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FIG. 15: Plot of average phonon occupation number for luq — 
2. Or, for several different bias conditions. Except for the most 
asymmetric case (fi^ = V s d), the phonon number diverges for 
large enough bias voltages. Also note that the onset voltage 
for the deviation of the phonon number from its equilibrium 
value (of zero) also continously shifts from /it = uo (most 
asymmetric bias) to \il = <^o/2 (symmetric bias) signaling 
inelastic cotunneling. 



Perturbative calculation for Phonon 
distribution function 



So far the out-of equilibrium Green's function tech- 
nique allowed us to calculate various averages such as 
current and mean phonon number. However it is also 
interesting to ask what the phonon probability distribu- 
tion function itself is under non-equilibrium steady state 
conditions. In order to do so we revert to the density ma- 
trix formalism developed in Section II. The complication 
in calculating this quantity is the nontrivial form of the 
operator equation Eg. 1191 One can however carry out the 
calculation for the rates to leading order in the electron- 
phonon coupling, using the diagrammatic language de- 
veloped by various authors2*2^. In implementing this we 
again find it convenient to be in the exact eigenstate ba- 
sis for the non- interacting system (Eq. I90|) . The leading 
order contribution to the rates is obtained by expanding 
the exponentials entering in the quantum rate equation 
H9l to leading order in the electron phonon coupling. The 
explicit form for the in-scattering and outscattering rates 
is therefore 

= Tr leads ^ dt{n\H e _ ph {Q)\n ± 1) (122) 
(n±l\H e _ ph (t)\n) 



He-ph 



k,k' ,a,b=L,R 



(123) 

On evaluating the above expressions we obtain the fol- 
lowing quantum rate equation for the distribution func- 
tion 



— — Pn' 1 (R n ->n+l + Rn->n-l) 
"T n _j_i J Vi+l— >n ~r * n _i -t^n— 1— m 

where the in-scattering rate is given by 

Rn^n+i =X 2 u; 2 (n+l)J £ 

E., fc L, ii r Jt/('-N)( 1 -/( e -^c»)) 

{(s-eo) 2 + ^}{( £ - £ o-^o) 2 + ^} 

and the out-scattering rate is given by 



(124) 



(125) 



= (n + l)r* 



Ea, b =i.,fl r ° r t./( e -M.)(l-/(e+wci-W,)) 
{(£-eo) 2 + ^}{(e-eo+^o) 2 + ^} 

We have used short hand notations r 



(126) 



nr out 



for the 



in I out 

cumbersome integrals that appear in the definition of 
etc. In terms of r in / out the quantum Boltzmann 
eciuation ll24l takes the form, 

P rh {{n+1) l™ +n} = pPl l{n+l)+{ n) P P^^L (127) 

Tout Tout 

It is easy to check that the above equation has a simple 
solution given by 



P? = (1 



Tout 



Tout 



(128) 



The quantity ■^ i2L - has been plotted for various combina- 
tions of r and cuo in Fig. El The figure illustrates that 
the rapidity with which the phonon distribution diverges 
with bias depends on the relative sizes of T and ujq. The 
stronger the coupling to the leads, the more easily the 
phonons equilibrate. 

In the high-T classical regime the rate equations for 
phonons for weak electron-phonon coupling has the same 
structure as Eq. 11271 but with modified scattering rates 
t 'in /out- (This has been explicitly shown in Appendix B). 
It therefore follows that even in the high-T regime, the 
phonon distribution function is given by Eq. 11281 

Note that on comparing Eqns. 112611251 with 
Ecins. [TD7I108I we find, not surprisingly, 



r m _^ H K - (U R - IT 4 ) 

r out n K + n R -u A 



(129) 



where the polarization LI were evaluated within the 
Keldysh Green's function approach. Thus in particular 
(N ph ) = = J2n nP n h , where D K has been 

calculated for phonons with no life-time broadening. 
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FIG. 16: Zero temperature conductance and current for 
equilibrated, unequilibrated and non-interacting levels under 
asymmetric bias conditions ({ir = 0,hl = V s d)- The phonon 
frequency uiq = 2.0r and r_L = Tjj = 0.5r. Moreover the 
electron-phonon dimensionless coupling strength is A = 0.25. 
Note that I (upper panel) is in units of and G (lower panel) 

is in units of 




^L /C0 

FIG. 17: Plot of ratio of phonon inscattering and outscat- 
tering rates for several different ratios of luq/T. The bias 
conditions are /i^ = — /j,r, and the onset of a non-zero r in at 
fiL = wo/2 is the sign of inelastic cotunneling. The high-T 
limit is approached for F <C <^o, when the onset of nonzero 
rin shifts to the sequential tunneling limit of fiL = u)q. 



V. LOW T TO HIGH-T CROSSOVER 

We would now like to connect our low-T quantum cal- 
culation with the high-T calculation discussed in Section 
III. The high-T limit may be reached by taking T > T 
in the Keldysh calculation, and the crossover from the 
T = to T ^> r case has been illustrated in figures IT81 
and 1191 for equilibrated and unequilibrated phonons re- 
spectively. The results for the equilibrated phonons have 
been presented for a rather large electron-phonon cou- 
pling Xcoq — 10r which strictly speaking is beyond the 
limits of validity of the perturbative approximation used 
here, in order to illustrate how the phonon sidebands 
evolve with temperature. 

The top panel of Fig. H8I illustrates how the elastic res- 
onance broaden with temperature, while the lower panel 
shows the broadening of the phonon side-band. As is 
evident from the two panels, at high temperatures the 
agreement between the high-T rate equation calculation 
and the quantum calculation is much better for the case 
of phonons with no life-time broadening. The effect of the 
life-time broadening due to interactions with electrons is 
to round off the phonon side-band further. 

Fig. ^5] illustrates the cross-over from the low-T 
to high-T regime under conditions of unequilibrated 
phonons and within the perturbative limit of Xluq = 0.5T 
and asymmetric bias. The phonon side-bands vanish for 
T > Xujq. 

VI. CONCLUSIONS 

A. Summary 

In this paper we have studied a simple model of an 
electron-phonon coupled quantum dot, involving a (pos- 
sibly degenerate) electron level coupled to leads and 
to phonons. The problem has four different important 
"internal" parameters: a dimensionless electron-phonon 
coupling A (defined for example in Eq.|2J, the ratio of the 
phonon frequency too to the broadening T of the on-site 
level due to coupling to the leads, the ratio of the tem- 
perature T to the level width T. (A fourth parameter, 
the ratio of the rate 7 eg at which the phonons relax to the 
heat bath characteristic of the device, to the mean elec- 
tron current flow rate has only been studied in limiting 
cases). 

The model admits two important sub-cases: of 
phonons coupled to the number of electrons on the 
molecule (Eq. [2J, and of phonons coupled to the dot- 
lead hybridization (Eq. [2J . Our formulation applies to 
both cases, but we have focussed mostly on the for- 
mer (McCarthy et al. have considered special features 
of the latter case in the classical regime). In this pa- 
per we have attempted to present a general framework, 
within which different special cases can be analyzed as 
desired. There are two important crossovers: the elec- 
tronic quantal-classical crossover controlled by T/Y, and 



18 




0.1 0.2 0.3 0.4 



> 




FIG. 18: Quantal-classical crossover effects in differential 
conductance spectra for equilibrated phonons. Upper panel: 
temperature dependence of zero bias "resonance" peak com- 
puted as described in text and compared to results obtained 
from classical rate equations (section III, but U = 0). Lower 
panel: temperature dependence of first phonon sideband, 
computed as described in text and compared to results ob- 
tained from classical rate equations. The rapid thermal 
smearing of both central peak and phonon sideband is evi- 
dent. Note that in the phonon side-band case, broadening 
of the phonon level due to electron phonon coupling leads to 
additional smearing not included in the rate equation model. 
The parameters are w = 100r, A = 0.1 and T = 1,5, 10r. 
The bias conditions are fiL = ~^r- Rate equation calcula- 
tion: A = 0.1, T = 10r. Note: dl/dV is in units of ^ZT- 



the phonon adiabatic/anti-adiabatic (phonon frequency 
long or short relative to inverse electron dwell time on 
molecule) crossover controlled by ujq/Y. (The limit of 
u>o <C T is not interesting) . In the classical limit (roughly, 
T/T > 1) this program has been carried out completely 
for all values of coupling A by us and by other work- 
ers. The relation between our results and those of other 
workers is discussed in detail in subsection B below. 

The model has two external "control" parameters; the 
source-drain voltage difference V sc i = hl—^r (see Fig. 1) 
and the molecule one-electron addition energy e' mea- 
sured relative to the average of the source-drain voltage 
(A'l + and also referred to as the gate voltage. It 

is well known that the conductance exhibits steps when- 
ever one of the lead chemical potentials passes through 
the level energy e' . The important new consequence of 
electron-phonon coupling is the appearance of steps such 
as those shown in Fig. [31 in the IV characteristics, when 
the source-drain voltage passes though an integer multi- 
ple of the phonon frequency. The existence of the phonon 
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FIG. 19: Quantal to classical crossover effects on IV curve 
for unequilibrated phonons. The parameters are u>o = 10r, 
A = 0.05 and the temperature for the classical calculation cor- 
responds to T = 10. Or. The quantum calculation has been 
performed for T = 0.1, 1, 10F. The results are for asymmetric 
bias fiR — 0.0, hl = V s d when the perturbation approxima- 
tion is valid. Note: I is in units of 

sidebands has been noted by several authors^. An im- 
portant new result (Fig. 0} of our work is the system- 
atic study of the dependence of the phonon heights on 
whether the phonon distribution function is controlled 
by the nonequilibrium current or is equilibrated to a heat 
bath (Aji et al. presented similar information for the co- 
tunneling regime). We have also studied (Figs. 1516115(1 
the non-equilibrium phonon distribution for different bias 
voltages. A further new result of our work is the theory 
for noise in the classical limit (section III C) showing that 
noise is a powerful spectroscopy of the degree to which 
the phonons are equilibrated (Fig. EJEJ, especially in the 
weak coupling limit. 

In the quantal limit (T < T), our treatment is re- 
stricted to low orders of perturbation theory in the cou- 
pling constant. Within this approximation we are able to 
obtain results (Figs El El EE! for the TV characteristics 
including both the "direct tunneling" (/i£ > e' > /ir or 
conversely) and "cotunneling" regimes {^l,^r > e' or 
conversely) and were able to treat the quantal to classi- 
cal crossover (Figs El El E3 1 ■ We presented both a di- 
agrammatic (Keldysh) calculation and a solution based 
on the construction of exact eigenstates, and confirmed 
that the peculiar broadening of the phonon distribution 
found in the "unequilibrated classical" case survives also 
in the quantum limit. 



B. Other work 

In this section we consider the relation of our results 
to the extensive existing theoretical literature. The sub- 
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ject was pioneered by Glazman and Shekter— who showed 
that the problem of a single electron transiting a res- 
onant level and coupled to phonons can be solved in 
complete analytic detail. They determined the form 
of the phonon side-bands in the transmission amplitude 
and showed how the resonant behaviour was modified 
by the electron-phonon interaction. Very similar results 
were subsequently obtained by Wingreen, Jacobsen and 
Wilkins^. However the single electron approximation 
used in these papers is not applicable to the case of leads 
containing a Fermi sea of electrons. The presence of other 
electrons blocks some of the intermediate states in the 
electron- phonon scattering process, changing the form 
of the eigenstates and crucially blocks some of the final 
states in the transmission process. This blocking ensures 
that the T = linear response conductance at resonance 
takes the quantized value if 2e 2 /h, whereas the on reso- 
nance transmission probabilities calculated in refsi^ are 
less than unity. This issue was very recently also dis- 
cussed by Flensbergi^. 

One crucial consequence of the presence of Fermi seas 
in the leads is the "floating" of the phonon side-bands in 
the electron spectral function and thus in the IV curve. 
At V s d = 0, the features in the spectral function occur 
at energies displaced from the fermi-level by integer mul- 
tiples of the phonon frequency thus the corresponding 
steps in the IV curve are observed when V s d is swept, 
but are not observed if V g (i.e., the mean lead Fermi 
level) is changed at fixed small V s d (see Fig. ITT|l . Sev- 
eral authors&iSiii employ approximations to the electron 
Green's function corresponding roughly to those of Glaz- 
man and Shekhter or Wingreen and coworkers. These 
miss the above physics and erroneously predicts phonon 
side-bands as V g is varied&iil. 

As also noted by Flensberg the approximations em- 
ployed by refsi^iiS amount to writing the dot Green's 
function as Gf(i) = e^ ie °- r ^{X^ (t)X(0)) (where the op- 
erators X have been defined in Eq. I31|) . This approx- 
imation becomes exact in the high temperature limit, 
and can be used as the starting point of an alternative 
derivation of the rate equations. However some of the 
literature^ who have used this approach have treated 
the X-operator matrix elements in an approximate man- 
ner which does not capture the structure giving rise to 
the step height variations displayed in our Fig. 

Finally Gogolin and Komnik have used an adiabatic 
(slow phonons, fast electrons) semiclassical approxima- 
tion to explore the strongly electron-phonon coupled 
regime. It is well known that at strong coupling the en- 
ergy as a function of phonon coordinate become bistable, 
signaling the onset of polaronic instability. Refii^ con- 
sidered the behaviour of the polaronic state under non- 
equilibrium conditions, and observed that a bistable I-V 
characteristic could result. By contrast our rate equa- 
tion analysis always leads to single valued I-V curves. 
The bistability corresponds to a first-order "energy land- 
scape" ; however the system under study is zero dimen- 
sional, and (at least within the approach used here) ther- 



mal and quantal fluctuations allow the system to explore 
all of phase space wiping out any phase transition be- 
haviour. The calculation of Gogolin and Komnik are 
presented for a different regime (T -C T), but we sus- 
pect that fluctuations would also eliminate the apparent 
transition in that case. An enhanced low frequency noise 
might however occur. 

Four recent paper have appeared which present results 
consistent with those presented here, but emphasizing 
somewhat different aspects of the physics. Mc Carthy 
et aliii used the rate equation approach of our section 
III, our results reproduce theirs; however the main inter- 
est if ref — was in the IV curves of phonons coupled to 
the dot-lead hybridization and the principal focus was on 
thermally equilibrated phonons. Our result extend theirs 
by treating the nonequilibrium case (see eg. Fig. |SJ and 
the noise (see eg. Fig. 01 ■ 

Fedorets and coworkers^ have analyzed the same 
problem as McCarthy et al. but in the limit of very 
weak coupling of phonons to a heat bath. They find 
that an instability occurs for V s d greater than a critical 
value. Interestingly in the weak coupling limit the critical 
V s d reported is identical to the critical V s d found in our 
calculation above which the phonon distribution become 
broad (see Fig. 0. According to the refsii& the instabil- 
ity occurs if one has both lead coupled and dot coupled 
phonons, whereas the broad distribution is generic. How- 
ever, this relation deserves further exploration. Here we 
observe that the calculations of refs^ (as well as those of 
our Section III) are based on a sequential-tunneling ap- 
proximation. In this approximation when the threshold is 
exceeded the distribution we find changes rapidly (in the 
weak electron-phonon coupling limit) from very narrow 
to very broad. However if cotunneling processes (in par- 
ticular the electron contribution to the phonon lifetime) 
are included, the transition becomes broadened with the 
onset moving to a lower voltage (see Fig. I15f) . We find 
the breadth of the distribution to depend strongly on the 
bias asymmetry , i.e., to the average occupation of the 
dot. 

Flensbergi^ has used an approximate equation of mo- 
tion approach to analyze the quantum limit, determin- 
ing in particular the equilibrium dot spectral function 
and presenting a clear analysis of the terms omitted in 
previous approaches^. Where there is overlap our re- 
sults are in agreement with his, however non-equilibrium 
conditions have not been considered in this reference^. 

Finally Aji, Moore and Varma — have considered 
phonon effects on the cotunneling spectrum. Their re- 
sults are in essential agreement with ours, however we 
note that the steps in dl/dV which they report to be 
vertical at T = are in fact smoothed by the phonon life- 
time arising from electron-phonon coupling (see Fig. I14|) 
or to coupling to a heat bath. Also we find that the 
steps are thermally broadened into indetectability at a 
relatively low T fixed by the step height and the phonon 
lifetime induced slope, rather than generically visible up 
to temperatures of order the phonon frequency, as stated 
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in the reference^. 

As far as experiments are concerned, three recent pa- 
pers^A have observed signatures in the current- voltage 
characteristics which are attributed to coupling to a 
molecular vibrational mode. A direct comparison to our 
results is not yet feasible because noise measurements 
have not been performed, and most of the samples stud- 
ied show one or at most two phonon sidebands. (Only 
one of the samples of Ref£ was reported to show more). 



C. Future Directions 

Finally we briefly mention a few issues raised by this 
work, but not resolved. One important research area is 
the extension of our quantum limit results beyond the 
perturbative regime (in particular to the nonequilibrium 
polaron limit), a second is to obtain the frequency de- 
pendence of the noise spectra in the classical and quan- 
tal regimes, a third issue is the crossover between equi- 
librated and uncquilibratcd phonons, and a fourth issue 
is to explore the connection between the nonequilibrium 
distributions we find and the instabilities discussed by 
the Chalmers groupie. Another interesting question is 
to explore the effect of image charges induced in the sur- 
rounding electrodes on I-V 26 which have been argued 
to be important in certain experiments-. Work in these 
directions is in progress. 
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APPENDIX B: HIGH-T LIMIT OF PHONON 
RATE EQUATIONS 



In this section we shall show how for weak electron 
phonon coupling and at high-T the quantum Boltzmann 
equation for phonons (Eq. 11241 11251 112611 reduces to the 
classical rate equations derived in Section II. For T ^> T, 
we can replace the lorentzian broadenings in Eg n il 2 5l and 
11261 by 5 functions i.e. ( ^s^jSyg ~~ * 7r5(a;)). Moreover 
for weak electron-phonon couplings we may identify the 
probability of single occupancy of the dot to be given by 
its value in the absence of electron phonon coupling, i.e. 
so that 
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Now the classical rate equations in Eq. I37l rewritten for 
the case of a single resonant level with U — take the 
form, 

P„° = -J* J2 R n°ra + £ J* (B3) 
a, 771 a,m 



APPENDIX A: ZERO TEMPERATURE PHONON 
POLARIZATION 

At zero temperature the function 1 — 2f(x) = sgn{x), 
and the integral Eq. 11091 and 11101 may be performed 
explicitly. The corresponding expressions are 
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where R^'^m e t c are defined in Eqns. [fH] and with 
U = 0. In the weak electron-phonon coupling limit, only 
transitions that change the phonon number at most by 
1 are allowed. In this limit the n+1 — » (n + l)A 2 T a . 
Moreover we may again factorize the joint probability 
distribution of having 0/1 electrons and m phonons as 
P^ 1 = P 0,1 P,p' 1 , After making these approximations in 
Eq. IB3I and IB4I and adding the two equations we obtain 



Pf = -P^iR^^+R^^) (B5) 

r pph p , pph p 



with the rates given by Eq. IB1I and IB 21 
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